Electrofreezing of Liquid Ammonia

Here we prove that, in addition to temperature and pressure, another important thermodynamic variable permits the exploration of the phase diagram of ammonia: the electric field. By means of (path integral) ab initio molecular dynamics simulations, we predict that, upon applying intense electric fields on ammonia, the electrofreezing phenomenon occurs, leading the liquid toward a novel ferroelectric solid phase. This study proves that electric fields can generally be exploited as the access key to otherwise-unreachable regions in phase diagrams, unveiling the existence of new condensed-phase structures. Furthermore, the reported findings have manifold practical implications, from the safe storage and transportation of ammonia to the understanding of the solid structures this compound forms in planetary contexts.

respectively. In addition, tests with a sample composed of an odd number of molecules (i.e., 101 (404 atoms)) were performed at a density of 0.80 g·cm −3 . The results reported in the main text refer to the simulation box containing 256 molecules at a density of 0.73 g·cm −3 .
As usual, in order to minimize spurious and nonphysical surface effects, the structures were replicated in space by employing periodic boundary conditions. Intensities of the electric fields were gradually increased with a step increment of 0.05 V/Å from zero up to a maximum of 0.70 V/Å. Whilst dynamics of 20 ps were performed in all zero-field cases, for each other value of the field intensity dynamics of 20 ps were ran for the simulated samples up to 128 molecules. For the biggest samples containing 256 NH 3 molecules, trajectories of 10 ps for each value of the field strength were accumulated. Tests aimed at checking the convergence of the presented results with longer simulations at a given field intensity (i.e., 100 ps) were executed. Besides, additional long simulations (≥ 100 ps) were ran to check the hysteresis upon switching off of the field. In all cases, solid-like structures returned into the liquid state within tenths of ps.

Wavefunctions of the atomic species have been expanded in TZVP basis sets with
Goedecker-Teter-Hutter pseudopotentials using the GPW method [3]. A plane-wave cutoff of 400 Ry has been imposed. Exchange and correlation (XC) effects were treated with the gradient-corrected Becke-Lee-Yang-Parr (BLYP) [4] density functional. In order to take into account dispersion interactions, the dispersion-corrected version of BLYP (i.e., BLYP+D3(BJ)) [5] was employed. Simulations were carried out at the nominal temperature of 252 K, similarly to other computational investigations [6]. In the AIMD simulations, the dynamics of the nuclei was simulated classically using the Verlet algorithm whereas the canonical sampling has been executed by employing a canonical-sampling-through-velocityrescaling thermostat [7] set with a time constant equal to 20 fs.

B. Path integral ab initio molecular dynamics
Path integral ab initio molecular dynamics (PI-AIMD) [2] simulations with a colored noise thermostat [8] of bulk liquid ammonia within a constant number, volume, and temperature (NVT) canonical ensemble were also performed with the aim of evaluating the role carried by Nuclear Quantum Effects (NQEs) in eventually preventing the electrofreezing observed in the classical nuclei AIMD simulations. NQEs may in principle play an important role in lowering the free-energy barrier for molecular dissociation and proton transfer [9]. The PI-AIMD sampling was obtained in the framework of the generalized Langevin equation thermostat (PIGLET) method [8] using six beads per atom. Such a thermostating technique, by accelerating the convergence, reduces the number of beads needed for mapping the classical system onto the quantum one [10]. All results refer to the centroid of the ring polymer beads. Being computationally more demanding than AIMD simulations, only the systems composed of 64 NH 3 molecules (i.e., 256 atoms) were simulated at the different densities reported in the previous section. Trajectories of 10 ps for each field intensity were gathered whereas 20-ps-long trajectories were accumulated in the zero-field regimes leading to a cumulative simulation time equal to ∼ 500 ps. An imaginary time-step of 0.5 fs for each bead of the PI-AIMD simulations was employed. Moreover, in the zero-field case, also tests employing imaginary time-steps equal to 0.4 fs and 0.3 fs were carried out.

C. Electric field implementation
The implementation of an external field in computations based on Density Functional Theory (DFT) can be achieved by exploiting the Modern Theory of Polarization and Berry's phases [11]. In a nutshell, the electrostatic interactions are determined by first-principles, as reported by Umari and Pasquarello [12]. For instance, they demonstrated that the functional where E KS is the Kohn-Sham energy functional, ϵ is the field intensity, and P is the polarization, is exploitable as energy functional for a variational approach to the finite-field problem as well. Fig. 1 of the main text and Fig. S1 here have been determined by means of the software TRAVIS [13,14] from the centres of the Maximally Localized Wannier Functions (MLWFs) [15,16] calculated on the fly during the ab initio molecular dynamics (AIMD) simulations. Molecular dipoles from MLWFs centers can be determined as:

Infrared (IR) spectra shown in
where e is the electron charge, r i is the position vector of the MLWF center i, Z j is the atomic number of the nuclei j whilst R j is the position vector of this latter. This way, the IR spectra at all the investigated field intensities of liquid ammonia were computed as the Fourier transform of the molecular dipole autocorrelation function along the simulation trajectories.
Experiments locate the IR vibrational bands ascribed to the symmetric and antisymmetric bending modes at 1050 and 1630 cm −1 , respectively [17]. As shown in Fig. 1 of the main text (black curve) and in Fig. S1 here, the employed Density Functional Theory (DFT) framework (i.e., BLYP+D3 [4,5]) localizes faithfully the center frequencies of the bending modes of liquid ammonia in the zero-field regime. In fact, whereas the frequency peak of the symmetric bending mode falls at 1080 cm −1 , the antisymmetric one is located around 1640 cm −1 , giving rise to a maximum error smaller than 3% with respect to the experimental data. Furthermore, the excellent agreement of our computations at zero-field regime with laboratory experiments results is also witnessed by the comparison of the location of the frequency peaks ascribed to the antisymmetric and symmetric N-H stretching bands.
In fact, whilst IR experiments of liquid ammonia locate the antisymmetric and symmetric stretching at 3380 and 3240 cm −1 [17], respectively, our first-principles simulations find peak frequencies equal to 3395 and 3255 cm −1 in this spectral region. This finding indicates that the treatment of the molecular dipoles correlations and the overall simulation setup (i.e., box size, trajectory length, etc.) -producing a discrepancy lower than ∼ 1% when compared to experiments -adequately describe subtle intra and intermolecular interactions in liquid ammonia. Thus, although neglecting Nuclear Quantum Effects (NQEs) in classical nuclei AIMD simulations generally results in blue-shifted stretching vibrations [18], also by virtue of errors compensation effects the BLYP+D3 functional appears to work particularly well for liquid ammonia.
As visible from the IR spectra of Fig. 1 of the main text and Fig. S1 here, the progressive growth of a high-frequency side mode of the low-frequency band is the manifestation of the additional topological constraint imposed by the field on the molecular dipoles. This vibrational Stark effect has recently been investigated in low-temperature (i.e., 10 K) ammonia molecules freely rotating in a Ar matrix and exposed to a DC electric field [19]. This way, ammonia molecules start librating rather than rotating under sufficiently intense electric fields. Besides, a contraction of the entire frequency range in the IR spectrum is recorded upon raising the field, an aspect appearing to characterize the general vibrational response of H-bonded liquids to externally applied electric fields. In other words, whereas the highfrequency bands are shifted by the field at lower frequencies, the low-frequency bands are blue-shifted, as shown in Fig. 1 of the main text and in Fig. S1.
Owing to the coupling between the external electrostatic potential gradient and the lo- value close to that recorded in bulk liquid water under static electric fields [21].
To track, at the molecular scale, the dynamical response of liquid ammonia to external static and homogeneous electric fields, the translational self-diffusion coefficient D of the ammonia molecules has been calculated from the mean square displacement (MSD) as where the indices i and j run on all pairs of first-neighbor molecules which at t 0 were H-  Van Hove correlation functions G(r, t) are particularly useful because their two-dimensional Fourier transform corresponds to the dynamic structure factor S(q, ω), which can be measured by neutron scattering experiments. To monitor molecular correlations we have calculated the temporal decay of the G 1 N N (t), a quantity measuring the correlations between first-neighbors nitrogen atoms over time (Fig. S5-c). Considering that when no correlations are present G 1 N N (t) ∼ 1, it follows that in the zero-field regime structural correlations between first-neighboring ammonia molecules vanish more rapidly than in presence of a finite field, as shown in Fig. S5-c. In addition, the field-induced enhancement of the spatial and temporal correlations is not merely referred to the nearest-neighbor species but applies at all scales, as visible from a direct comparison of the Van Hove correlation function between nitrogen atoms G N N (r, t) at zero field ( Fig. S5-a) and upon applying a field intensity of 0.50 V/Å (Fig. S5-b). Typical intermolecular correlations occurring at different distances persist, indeed, over the simulated timescales in presence of the field, as visible from the red/yellow/orange regions in the contour map of the Van Hove correlation function plotted in Fig. S5-b. By contrast, a diffuse green area associated with G N N (r, t) ∼ 1 testifies the rapid temporal decay of the spatial structural correlations in liquid ammonia in absence of the field due to the presence of weak H-bonds (Fig. S5-a).
The Voronoi cell of individual ammonia molecules was widely used to characterize local structures in liquid systems. Asphericity of the Voronoi cell parameter η features the shape of the resulting Voronoi cell and is given by the following expression: where V is the Voronoi cell volume and A is its surface area. Asphericity is here calculated based on the Voronoi analysis on the nitrogen atoms of ammonia. This unit-less parameter is proving that the field acts as an order-maker agent in the ammonia structure (see Fig. S11).
On the other hand, it is worth remarking that the orientational order parameter q 6 is far more sensitive to the field-induced structural changes than the asphericity of the Voronoi To quantitatively evaluate growing spatial correlations in the simulated systems we hence employed the translational order parameter t defined as [23]: In this expression density-density correlations are detected by integrating over the absolute value of the total correlation function h(ξ) = g(ξ) − 1, where g(ξ) is the RDF determined as a function of the rescaled spatial coordinate ξ = rρ 1/3 and ξ c is a numerical cutoff imposed by the finite size of the simulation box (i.e., ξ c ∼ 3 in our case). The strength of this choice resides in the fact that it represents a measure of the order which is independent from the specific crystalline symmetry, this latter being a priori unknown. Fig. S10-b shows the trend of t as a function of the field strength.
Aside from translational order, crystalline structures are featured by orientational order. Let us consider the global orientational order parameter defined as [24]: where Q lm = ⟨Y l,m (θ i,j , ϕ i,j )⟩ is the average of the spherical harmonics Y l,m of all bonds connecting the nitrogen atoms in the system. q 6 is the orientational order parameter obtained by imposing l = 6 and it was chosen in the current investigation to monitor the field-induced increase of the orientation order in liquid ammonia (see Fig. S10-a). hydrogen atoms appear to be disordered.